


library(mvtnorm)
library(arm)
library(BRugs)
library(R2OpenBUGS)
library(coda)
library(car)
library(foreign)
library(vioplot)

library(RItools)
library(xtable)

# Change this to your working directory.  Note this path makes Dropbox work at least on a windows machine
setwd("~/../Dropbox/iREDS/Data/Analysis/single_items")


dataset<-read.dta(file.choose()) # read in a Stata file "ireds_dataset_wide_v12.dta"
# the next line is if you want to replicate the complete case MLE models -- we don't use!
####dataset1<-dataset[complete.cases(cbind(dataset$t0_genrcr_discusschangedviews, dataset$t1_genrcr_discusschangedviews)),]


xBalance(condition~urm+male+pi_role , data=dataset, report=c("all"))




# Data preprossing for the Congdon model 


# create the matrix of lab adjacencies

LabAdj=matrix(data=0, nrow=length(dataset$labid), ncol=length(dataset$labid))
for (i in 1:length(dataset$labid)) {
  for (j in 1:length(dataset$labid)) {
    if (i!=j) { 
      if (dataset$labid[i]==dataset$labid[j]) {
        LabAdj[i,j]<-1
      }}}}



#### Create Vector of labmate mappings to read into Congdon programs: 

LabMap<-list(map=NULL)
for (i in 1:length(dataset$labid)) {
  for (j in 1:length(dataset$labid)) {
    if (LabAdj[i,j]==1) LabMap<-list(map=c(LabMap$map, j))
  }}

LabC<-list(C=0)
nadj<-0
for (i in 1:length(dataset$labid)) {
  nadj<-sum(LabAdj[i,])+nadj
  LabC<-list(C=c(LabC$C, nadj))
}
map <- LabMap$map
C <- LabC$C

one<-rep(1, length(dataset$labid))
numLabNeigh<-sum(as.numeric(t(LabAdj%*%one)))
Cn <- numLabNeigh





## create data 


# set constants 
n_respondents <- length(dataset$labid)
n_responses <- 4
n_labs <- max(dataset$labid) # not used in model 
labid <- dataset$labid # not used in model 


 

# create matrices of the pre and post outcomes.  You can include an arbitrary number of items/questions provided
# you have at least three to identify the latent variables.
S1.0 <- as.matrix(cbind(dataset$t0_genrcr_discusschangedviews))
S2.0 <- as.matrix(cbind(dataset$t0_authorpolicy_yesno, dataset$t0_datamanagement_yesno))

S1.1 <- as.matrix(cbind(dataset$t1_genrcr_discusschangedviews))
S2.1 <- as.matrix(cbind(dataset$t1_authorpolicy_yesno, dataset$t1_datamanagement_yesno))

S3 <- as.matrix(cbind(dataset$m1_dmsystem,	dataset$m1_ethicstraining,	
                      dataset$m2_dmsystem,	dataset$m2_ethicstraining))

condition <- as.numeric(dataset$condition)

cuts<-c(-1,0,1)



data<-list("S1.0", "S2.0", "S1.1", "S2.1", "S3",
           "condition", "n_respondents", "n_responses", "C", "map", "Cn") 

inits <- function() {
  list(S1.beta.1=rnorm(1), S1.beta.2=rnorm(1), S2.beta.1=rnorm(2), S2.beta.2=rnorm(2), 
       S1.beta.0=c(cuts), S1.beta.00=c(cuts),  S2.beta.0=rnorm(2), S2.beta.00=rnorm(2),
       S3.beta.2=rnorm(4), S3.beta.0=rnorm(4),
       S1.rho=0.1, S2.rho=c(0.1, 0.1), S3.rho=c(0.1, 0.1, 0.1, 0.1),
       S1.e=rnorm(n_respondents), S2.e=rmvnorm(n_respondents, mean=rep(0, 2)), S3.e=rmvnorm(n_respondents, mean=rep(0,4))
  )} 


# Output the data and inits files to the working directory for the analysis...
setwd("~/../Dropbox/iREDS/Data/Analysis/single_items")
bugs.data(data, dir=getwd(), digits=5, data.file="data.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits1.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits2.txt")
bugs.data(inits(), dir=getwd(), digits=3, data.file="inits3.txt")



# Read in OpenBUGS results and create figure


coda.data<-read.openbugs(stem="")
results<-as.data.frame(rbind(coda.data[[1]], coda.data[[2]], coda.data[[3]]))


nrows<-length(results[,1])
ncols<-3 # this is the number of rows in the figure
figdata<-as.data.frame(matrix(rep(NA, nrows*ncols*2), nrow = nrows*ncols, ncol=2))

for (j in 1:ncols) {
  for (i in 1:nrows) {
    if (j==3) {
      figdata[i+2*nrows,1] <- results$S1.beta.2[i]
      figdata[i+2*nrows,2] <- j}
    if (j==2) {
      figdata[i+nrows,1] <- results$S2.beta.2.2[i]
      figdata[i+nrows,2] <- j}
    if (j==1) {
      figdata[i,1] <- results$S2.beta.2.1[i]
      figdata[i,2] <- j}
    }
}

figdata[,2] <- factor(
  figdata[,2],
  levels = c("1", "2", "3"),
  labels = c("Lab Has Authorship Policy", "Lab Has Data Management Plan", "Discussion Changed Ethics Views")
)




library(ggplot2)
library(ggridges)
library(ggpubr)
library(viridis)

theme_set(theme_ridges())



## USE THIS ONE -- THIS IS THE VERSION IN THE PAPER!
ggplot(figdata, aes(x=figdata[,1], y=figdata[,2], fill=0.5 - abs(0.5-..ecdf..))) +
  stat_density_ridges(geom = "density_ridges_gradient", scale=0.9, calc_ecdf = TRUE) +
  scale_fill_viridis(name = "Prob(Tail)", direction = -1) +
  labs(
    title = 'Treatment Effect Estimates',
    subtitle = 'Log Odds Scale',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())





ggplot(figdata, aes(x = figdata[,1], y = figdata[,2], fill = ..x..)) +
  geom_density_ridges_gradient(scale = 0.9, rel_min_height = 0.01, gradient_lwd = 1.) +
  scale_x_continuous(expand = c(0.01, 0)) +
  scale_y_discrete(expand = c(0.01, 0)) +
  scale_fill_viridis(name = "Scale", option = "D", begin=0, end=1, direction = 1) +
  labs(
    title = 'Treatment Effect Estimates',
    subtitle = 'Log Odds Scale',
    x = "Posterior Distribution"
  ) +
  theme_ridges(font_size = 13, grid = TRUE) + theme(axis.title.y = element_blank())


op <- par(no.readonly = TRUE) # the whole list of settable par's.
par(mar=c(5,15,5,1), mfcol=c(1,1), xpd=FALSE, las=1)
# Draw the plot
with(results , vioplot(S1.beta.2, S2.beta.2.2, S2.beta.2.1, #S3.beta.2.1, S3.beta.2.2, S3.beta.2.3, S3.beta.2.4,  
                       col=topo.colors(3), range = 1.75, lwd=2, border = NA, lineCol = "gray",
                         names=c("Discussion Changed Ethics Views", "Lab Has Data Management Plan", 
                                 "Lab Has Authorship Policy"
                                  #"Learned Online Data Management - Midpoint 1", "Received Ethics Training - Midpoint 1",
                                 #"Leaned Online Data Management - Midpoint 2", "Received Ethics Training - Midpoint 2"
                                 ), 
                        main=paste("Treatment Effect Estimates (Log Odds Scale)"), xlab=paste("Posterior Distribution"),
                         horizontal = TRUE))  
abline(v=0)
#abline(h=3.5, lty=2)
#text(x=-1, y=5.7, "Lab Communication")
#text(x=-1, y=5.4, "Scales")
#text(x=-1, y=2.7, "Training Content")
#text(x=-1, y=2.4, "Scales")
par(op)



# results:
length(results$S1.beta.2[results$S1.beta.2>0])/length(results$S1.beta.2)
length(results$S2.beta.2.1[results$S2.beta.2.1>0])/length(results$S2.beta.2.1)
length(results$S2.beta.2.2[results$S2.beta.2.2>0])/length(results$S2.beta.2.2)





